Convergence of an s-wave calculation of the He ground state. 
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The Configuration Interaction (CI) method using a large Laguerre basis restricted to £ = orbitals 
is applied to the calculation of the He ground state. The maximum number of orbitals included was 
60. The numerical evidence suggests that the energy converges as AE N « A/N 7 ^ 2 + B/N s ^ 2 + . . . 
where /V is the number of Laguerre basis functions. The electron-electron 5-function expectation 
converges as A8 N w A/N 5 ^ 2 + B/N a ^ 2 + . . . and the variational limit for the 1 = basis is estimated 
as 0.1557637174(2) <jq. It was seen that extrapolation of the energy to the variational limit is 
dependent upon the basis dimension at which the exponent in the Laguerre basis was optimized. 
In effect, it may be best to choose a non-optimal exponent if one wishes to extrapolate to the 
variational limit. An investigation of the Natural Orbital asymptotics revealed the energy converged 
as AE N » A/N 6 + B/N 7 + . . . while the electron-electron S- function expectation converged as 
A5 N ~ A/N 4 + B/N 5 + . . .. The asymptotics of expectation values other than the energy showed 
fluctuations that depended on whether N was even or odd. 

PACS numbers: 31.10.+Z, 31.15.Pf, 31.25.Eb 
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I. INTRODUCTION 

There have been a number of studies of the conver- 
gence of the configuration interaction (CI) expansion of 
the helium ground state QBHESIlSlll followin g 
the pioneering work of Schwartz These studies have 
investigated the convergence of the energy with respect 
to the number of partial waves included in the wave func- 
tion and also with respect to the dimension of the radial 
basis. 

It has been known since 1962 Q that the energy con- 
verges slowly with respect to J, the maximum angular 
momentum of any orbital included in the CI expansion. 
In particular the leading term to the energy increment is 
known to behave as 



AE J = (E) J - {E)- 1 - 1 



At 
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at high J. Later work 0, SH H showed that the 
energy increments can be written more generally as 
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where explicit expressions for Ae and Be exist, namely 

A E = -6tt 2 J |#(r,r,0)|Vdr = -0.074226 (3) 

B E = J |*( r ,r,0)| 2 r 6 dr = -0.030989 . (4) 

No expressions for Ce exist. The numerical values in 
eqs. © and are obtained from close to exact wave 
functions 0- 

However, the convergence with respect to J represents 
only one aspect of the convergence problem. Just as im- 
portant is the convergence with respect to the dimension 
of the radial basis N, for a given J. How do the incre- 
ments to E with increasing N 
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behave? In effect, what are the values of p and i? This 
aspect of the CI expansion is not as well understood as 
the convergence with J and there have been no studies 
equivalent in sophistication to those of Schwartz [9] , Hill 
and Kutzelnigg and collaborators 0, @- Some at- 
tention has been given to the radial convergence of the 
hydrogen atom in gaussian basis sets [lOj. The seminal 
investigation of Carroll and collaborators concluded that 
p Ri 6 for a natural orbital (NO) basis 0,^3- This result 
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has been quite influential, and can be regarded as ulti- 
mately motivating the use of principal quantum number 
expansions to extrapolate energies to the infinite basis 
limit from correlation consistent basis sets El ■ More re- 
cently, Goldman performed a regression analysis to give 
p»5.7 for a NO basis and p«3.8 for a Slater basis with 
a common exponent [l3T |. 

The radial basis sets used for the configuration inter- 
action or many body perturbation theory treatments of 
atomic structure can be broadly divided into two classes. 
In the first class, one defines a box and a piece- wise poly- 
nomial (e.g a spline) is used to define the radial depen- 
dence of the wave function in the interior of the box. The 
properties of the radial basis are determined by the size 
of the box, the number of knot points, and where they are 
located. The other approach typically expands the wave 
function in terms of a basis of functions with a convenient 
analytic form, examples would be an evenly tempered set 
of Slater type orbitals (STOs) 0, Q (this type of basis 
set is often optimized with respect to a couple of param- 
eters used to defined a sequence of exponents) or a set of 
Laguerre type orbitals (LTOs) [H EE 13 

The two most recent examples of these two approaches 
arc the calculations by Decleva et al (B-splines) Q and 
Sims and Hagstrom (Slater basis) [7j which are the 
biggest calculations of their respective type. The B-splinc 
calculation has given estimates AE J increments that are 
believed to be accurate to within hartree or better. 
One of the reasons this accuracy is possible is that AE J 
varies smoothly as the number of knot points is adjusted. 
This made it possible to obtain reasonable estimates of 
the infinite basis limit. Their estimate of the s-wave limit 
was accurate to better than 10 -9 hartree. Achieving this 
extreme level of accuracy was not possible when using 
the Slater basis [?J since linear dependence issues made it 
problematic to expand their radial basis to completeness. 
Indeed, resorting to REAL*24 arithmetic still resulted in 
an error of 4 x 1CP 6 hartree. 

In some investigations of the convergence properties 
of the CI expansion for the helium atom Q and mixed 
electron-positron systems 0| ^ became apparent that a 
better understanding of the issues that influence the con- 
vergence of the radial basis was desirable. For example, 
it was apparent that the dimension of the radial basis 
should be increased as J increases in order to ensure the 
successive AE J increments are computed to the same 
relative accuracy |8l Il7l Il8j . In addition, it was readily 
apparent that extrapolation of the radial LTO basis to 
the N — ► oo limit was not straightforward. 

In this work, we investigate the radial convergence of 
the CI expansion for a more manageable model of the 
helium atom with the orbitals restricted to the £ = 
partial wave. The linear dependence issues that are such 
a problem for a Slater basis are eliminated by choosing 
the radial basis to consist of LTOs 

(formally, the 

LTO basis spans the same space as the common exponent 
Slater basis, i.e. r n * exp(— Ar)). We note in passing the 
previous work of Hol0ein [16| who also investigated the 



convergence of a (small) LTO basis for an s-wave model 
of helium. Initially, we examine the merits of using an 
LTO basis with the exponent optimized to the basis di- 
mension. The nature of the asymptotic expansion for the 
energy increments is then deduced. Finally, the density 
matrix for our best wave function is diagonalized and the 
convergence properties of the natural orbital expansion 
are also determined. As part of this analysis, attention 
is also given to the convergence of the electron-electron 
coalescence matrix element since it arises in calculations 
of the two-electron relativistic Darwin correction p(J and 
electron-positron annihilation pi) . 



II. THE s-WAVE ENERGY FOR HELIUM 

The non-relativistic hamiltonian for the 1 S e ground 
state of helium 
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is diagonalized in a basis consisting of anti-symmetric 
products of single electron orbitals 

= 0) =£cy A12 (* W */ii|OO)0i(ri)&(r a ). (7) 



The functions 4>(r) are single electron orbitals Laguerre 
functions with the radial form 

Xa(r) = N a r e exp(-A Q r)L^+ 2 l 1 (2A a r) , (8) 

where the normalization constant is 



N a = 



The function L^+|l 1 (2A a r) is an associated Laguerre 
polynomial that can be defined in terms of a confluent 
hypergeometric function E,H3, Hj|. In the present work 
£ is set to zero. The basis can be characterized by two 
parameters, N, the number of LTOs, and A the exponent 
characterizing the range of the LTOs. It is normal to use 
a common exponent, A and when this is done the basis 
functions form an orthogonal set. The exponent can be 
optimized to give the lowest energy and when this is done 
one can define Am to be the value of optimal A for a basis 
of dimension M. The value of Am was seen to increase 
M increased. 

Calculations for basis sets with Am optimized for di- 
mensions of 10, 20, 30, 40, 50 and 60 have been per- 
formed. The optimal exponents are listed in Table[l]along 
with the energy for those values of M. The calculations 
were expedited by using the result that the two-electron 
Slater integrals for any A are related to each other by 
a very simple scaling relation. Accordingly, the list of 
Slater integrals could be generated once (by numerical 
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TABLE I: Comparison of different CI calculations of the s-wave model of the He atom ground state. The expectation value of 
the electron-electron 5-function (in ag) is denoted as {5}. The data in the {E) M and (E) 60 columns are the energies (in hartree) 
with N = M and N = 60 basis sets respectively. The data in the {E°°) and (5)°° columns are obtained by doing an explicit 
calculations with N = 60 and then adding in the 60 — » oo correction assuming an A/N~ p asymptotic form. 



M Am 


(E) M 


(Ef° 


(E)°° 


(S) M 


(d)°° 


10 3.07 


-2.879 022 691 296 


-2.879 028 727 964 


-2.879 028 7667 


0.155 922 600 334 


0.155 763 879 


20 4.80 


-2.879 028 507 141 


-2.879 028 754 899 


-2.879 028 7671 


0.155 789 345 524 


0.155 763 796 


30 6.45 


-2.879 028 726 601 


-2.879 028 761 447 


-2.879 028 7672 


0.155 772 304 341 


0.155 763 771 


40 8.04 


-2.879 028 756 467 


-2.879 028 763 935 


-2.879 028 7673 


0.155 767 637 040 


0.155 763 760 


50 9.57 


-2.879 028 763 441 


-2.879 028 765 118 


-2.879 028 7670 


0.155 765 847 239 


0.155 763 757 


60 11.10 


-2.879 028 765 650 


-2.879 028 765 650 


-2.879 028 7661 


0.155 765 002 122 


0.155 763 860 


Decleva el al [4] 




-2.879 028 767 289 








Goldman el al 




-2.879 028 767 319 









integration using gaussian quadrature 15]) for a given A, 
and then recycled by rescaling for calculations at differ- 
ent A (refer to Appendix B). Calculations with TV rang- 
ing from 1 to 60 have been performed for the six values 
of Am listed in Table [I] The dimension of the hamilto- 
nian for the largest calculation was 1830. The quanti- 
ties listed in the tables and the text are given in atomic 
units. The most precise energy for the helium s-wave 
model is that of Goldman [jjj, |24j who used a basis writ- 
ten in terms of r< , r> co-ordinates to obtain an energy 
ofE = -2.879028767319214 hartree. 



III. SIMPLE POWER LAW DECAY 



All observable quantities can be defined symbolically 



as 



JV 



(x) n = j2 Axn 



(10) 



where AX™ is the increment to the observable that oc- 
curs when the basis dimension is increased from n — 1 to 
n, e.g. 

AX" = (X) n - (X)"- 1 . (11) 
Hence, one can express the limiting value formally as 



A. Use of quadruple precision arithmetic 



The present calculations were all performed with 
quadruple precision arithmetic. It was only possible to 
get energies precise to 13 significant digits for the largest 
calculations when double precision arithmetic was used. 
This was caused by roundoff error gradually accumulat- 
ing during the course of the rather extensive calculations 
and the 13 digits appears to be the limit that can be 
achieved for double precision arithmetic (some experi- 
mentation revealed that the last 2 digits of the 15 double 
precision digits were sensitive to different Fortran com- 
pilers and even the optimization options of those com- 
pilers). The analysis requires investigation of the energy 
differences of eq. (JjJJ, and these energy differences can 
be rather small (e.g. AE 60 — 1.5 x 10~ 10 hartree for 
the Aeo basis). The fluctuations caused by roundoff did 
have a noticeable impact on the parameters derived from 
these energy differences at large N. These fluctuations 
were removed once quadruple precision arithmetic was 
adopted. 



(x) = (x) N + J2 Axn 



(12) 



n=N+l 



The first term on the right hand side will be determined 
by explicit computation while the second term will be 
estimated. Obtaining an estimate of the remainder term 
does require some qualitative knowledge of how the AX N 
terms decay with N. For example, previous computa- 
tional investigations indicate that the natural orbital de- 
composition leads to a AE N N~ 6 dependence [H |25|. 
As far as we know, there have not been any detailed in- 
vestigations of the N dependence of a Laguerre basis. 

A useful way to analyze the convergence is to assume 
the increments obey a power law decay of the form 



AX 



Ax 
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and then determine the value of p from two successive 
values of AX using 



P = In I — : „»r ) /I" 



V AX 



N 



N 



N - 1 



(14) 



Figure ^ plots the exponent derived from the energy 
increments for six different values of Am . The succession 
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FIG. 1: The exponents pe as a function of AT for the s-wave 
calculations of the He ground state. The different curves were 
obtained from the LTO basis sets with the exponents listed 
in Table □ 

of curves show that pe tends to peak at values larger than 
10 at an intermediate AT and then shows a tendency to 
decrease. The value at N — 60 was pe ~ 3.7 for the most 
of the curves shown in figure Q 

The salient point to be extracted from Figure ^ is that 
the value of pe for a given Am at N = M is quite dif- 
ferent from the asymptotic value, e.g. the value of pe 
for the A20 curve is much larger at N = 20 than it is 
at N = 60. This is quite an annoying result. Ideally, 
one would like to perform the largest calculation with 
the exponent optimized for that dimension basis. Then 
the specific form of the power law decay would be esti- 
mated by analyzing the energies obtained from a series 
of slightly smaller calculations. This information would 
subsequently be used to estimate the energy or other ex- 
pectation value in the variational limit. However, this 
is not possible since the energy increments will not have 
achieved their asymptotic form. 

Although there are problems in using an optimized ex- 
ponent, it may still be possible to analyze a sequence of 
energies from a calculation with a non-optimized expo- 
nent and thereby estimate the variational limit. Assum- 
ing that the increments obey eq. , one can write 

A x = N p AX N , (15) 

and thus the n > N remainder term 

can be derived from (X) N - 2 , (X)^ 1 and (X) N [l[l7|. 
When this remainder was evaluated in this work, the first 
10000 terms of the sum over n were computed explicitly. 
Then the approximate relation eq. (|16f) was used. 



FIG. 2: The extrapolated N — * 00 limit for the He ground 
state energy for different values of A a/. The exact s-wave 
energy is taken from the J = calculation of Goldman fl3jl . 

Figure |3 shows the estimated variational limit as a 
function of N for the A^ listed in Table |U An explicit 
calculation including A" LTOs was initially performed to 
determine (E) N . Then eq. I|16|) was used to estimate the 
remainder and hence deduce the variational limit. The 
variational limits in Tabled were extracted from the cal- 
culations with A^ = 60. The exact variational limit can 
be predicted to the 9th digit after the decimal point. The 
most inaccurate estimate of the variational limit is that 
from the Ago calculation. So the calculation that is ex- 
plicitly optimized at A^ = 60, (i.e. with and gives 
the best energy at A" = 60, gives the worst estimate of 
the variational limit! 

A CI calculation of the Li + ground state restricted 
to the I — partial wave was also performed to check 
whether the conclusions above were peculiar to He. Once 
again, the exponent pe was approximately 3.7 at A^ = 60 
and the convergence pattern for N < M was distorted 
by optimization of the exponent. 



IV. THE 5-FUNCTION EXPECTATION VALUE 

Part of the motivation for the present study is to gain 
a better understanding of how to perform CI calcula- 
tions for mixed electron-positron systems. Apart from 
the energy, the next most important expectation value 
for a positronic system is the electron-positron annihila- 
tion rate |2l|. The annihilation rate is proportional to 
the expectation of the electron-positron delta function, 
and has the inconvenient property that it is even more 
slowly convergent than the energy [TEl |26| . Accordingly, 
the convergence of the electron-electron (5-function is in- 
vestigated using the methodology previously used for the 
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Although distortions in the convergence pattern are 
still present, they are less severe than the energy since 
the successive AS N increments are larger. As a rule, ps 
was at least 10% larger than 2.6 at N = M . A choice 
of N > (M + 10) would generally lead to ps being in the 
asymptotic region. 

Figure 0] shows the estimated variational limit of (5)°° 
as a function of N for the six different values of Am 
listed in Table [I] An explicit calculation including N 
LTOs was initially performed, then eq. I)16|l was used 
to estimate the variational limit. A variational limit of 
(S) = 0.1557637174(2) Oq (see later discussion) was as- 
sumed for plotting purposes. The notable feature is that 
the Aeo estimate of the limit at N = 60 is one of the least 
accurate. 



V. A CLOSER LOOK AT THE ASYMPTOTIC 
POWER LAWS 



FIG. 3: The estimated exponents ps as a function of N for 
the LTO calculations of the He ground state (8) . The different 
curves were obtained with the LTO basis sets listed in Fable 
□ 



energy. The only independent investigation of this quan- 
tity for an s-wave model of helium was by Halkier et ol 
who obtained 0.155786 afj. 
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FIG. 4: Fhe convergence of (6) for the He ground state as 
a function of N. Fhe absolute value of the difference of the 
extrapolated N — > oo limit subtracted from 0.1557637174 Oq 
is plotted. 

Figure|3]tracks the behavior of the exponent ps derived 
from eq. lfl4*l) . It can be seen that ps achieves values ex- 
ceeding 10 before it decreases to its asymptotic values. 
The present calculations give ps ~ 2.6 at N = 60 al- 
though ps is still exhibiting a slow but steady decrease. 



Figures |S] and show the behavior of the pe and ps 
versus — j= for values of TV greater than 16. Curves are 
not shown for all the Am exponents. In some cases, the 
values of ps did not fall within the plotting window. 

The notable feature to be gleaned from both set of 
curves is the essentially linear behavior Pdeita with respect 
to — = for values of N greater than 16 and a visual in- 
spection suggests that the limiting exponents is ps = 2.5. 
The purely visual evidence that the lowest order term of 
Pe is 0(N^ 1 ^ 2 ) is not as compelling as that for Pdeita, 
but by analogy with this form has been assumed and 
subjected to extensive testing. 

More substantial evidence is provided by a ht of the p 
vs N data to an inverse power series of the form 



P = Pa 



N p 

£ 

i=l 



(17) 



What we have done is fit (N p + 1) successive pe or ps 
values to eq. i|17fl for the Aio data sequence. The results 
of those fits are given in Table ITT1 

Using a 4 or 5 term fits (and 6 term fit in the case 
of pe) results in limiting exponents very close to either 
3.5 (for pe) or 2.5 (for ps). These estimates were also 
reasonably stable. For example, the value of pe for a 
5-term fit for a data sequence terminating at N = 50 (as 
opposed to N = 60 in Table HTJ) was 3.4970. 

The validity of the series, eq. i|17|) immediately suggests 
that the asymptotic forms for AE N are 
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(In Appendix A it is demonstrated that an exponent vari- 
ation of p = po + B/y/N arises from an inverse power se- 
ries in AX N with a leading term of B/N Po and with the 



TABLE II: Results of using 2, 3, 4, 5 or 6 term inverse power 
series (corresponding to N p — 1,2,3,4 and 6) to determine 
the limiting value of the exponents and the energy and 5- 
function (in <Iq). The results for the LTO data sequences 
were taken from the largest calculations for the Aio exponents. 
The results for the NO data sequences were extracted at N = 
20 while those for the energy optimized LTO sequence were 
determined at N = 30. Data entries with an asterisk, *, were 
obtained using a weighted average (as described in the text) 
due to fluctuations depending on whether the N was even or 
odd. 





P0 


Pi 


A x 


(X)°° 






LTO data 


sequence: Aio 


basis 








(E) 




1 


3.4310 


2.0481 


-0.00182 


-2.879 028 767 0519 


2 


3.5121 


1.1301 


-0.00225 


-2.879 028 767 3496 


3 


3.4978 


1.1267 


-0.00212 


-2.879 028 767 3154 


4 


3.4988 


1.0958 


-0.00215 


-2.879 028 767 3196 


5 


3.5012 


1.0032 


-0.00215 - 


2.879 028 767 31920 








(5) 




1 


2.4809 


0.9444 


-0.00535 


0.155 763 7540 


2 


2.4975 


0.6872 


-0.00562 


0.155 763 7156 


3 


2.5024 


0.5760 


-0.00560 


0.155 763 7172 


4 


2.4989 


0.6810 


-0.00559 


0.155 763 7174 






NO data sequence 










(E) 




1 


5.9916 


2.7712 


-0.2959 


-2.879 028 767 3054 


2 


5.9971 


2.5553 


-0.2995 


-2.879 028 767 3176 


3 


5.9959 


2.6256 


-0.2996 


-2.879 028 767 3177 


4 


5.9671 


4.7513 


-0.2999 


-2.879 028 767 3179 








(5) 




1 


3.9973 


1.620* 


-0.1957 


0.155 763 7197 


2 


3.9980* 


1.681* 


-0.1962* 


0.155 763 7177* 


3 


3.9997* 


1.599* 


-0.1963* 


0.155 763 7175* 




Energy optimized LTO data sequence 








(E) 




1 


5.6562 


149.13 


-1.543 


-2.879 028 767 3333 








(S) 




1 


3.8093* 


0.3839* 


-0.3879* 


0.155 763 7191* 



power increasing by ^J~N for successive terms) . Although 
eqs. <|18fl and l|19|) are best described as a conjecture, the 
numerical evidence in support of the conjecture will be 
seen to be strong. 

The applicability and utility of eqs. (|T%|) and l(T§|l was 
tested by fitting these equations to (E) N and (5} N values 
and then using eq. 1161) to determine the N — > oo limits 
for the individual terms. Asymptotic series with up to 
5-terms (i.e. N p = 4) were also investigated. 

Pigure|7|shows (E)°° for the Aio basis using the asymp- 
totic series of different lengths to estimate the N — > oo 
correction. It is noticeable that all the representations 
of eq. (|18|) exhibit better convergence properties than 
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FIG. 5: The exponents pe as a function of for the s-wave 
calculations of the He ground state. The different curves were 
obtained from the LTO basis sets with the exponents listed 
in Table □ 

eq. <|13|) and the 6-term representation has the best con- 
vergence characteristics for A^ > 30. (It should be noted 
that the 3-, 4-, 5- and 6- term extrapolations to (E)°° ex- 
hibited fluctuations of order 10~ 5 to 10~ 9 hartree when 
the calculation was performed in double precision arith- 
metic). The increasingly better convergence character- 
istics as N p increases is consistent with eq. (|18f) be- 
ing the asymptotic form describing the energy conver- 
gence with respect to a LTO basis. The estimated (E)°° 
limits at A^ = 60 for the various asymptotic expres- 
sions are given in Table [H] The 6-term estimate was 
(E) 00 = -2.87902876731920 hartree which agrees with 
the value of Goldman, namely E = -2.87902876731921 
by 1 x 10~ 14 hartree. The precision of the extrapolated 
value exceeds the precision of raw (E) 60 energy by a fac- 
tor of 1,000,000! This improvement is best placed in per- 
spective by noting that the Aio calculation would have 
to be extended to iV w 10 6 to achieve the same level of 
precision. 

This extreme accuracy is not reproduced if one uses 
other forms for the asymptotic series. For example, mak- 
ing the choice AE N = A E /N 4 + B E /N 5 + ... results 
in much poorer estimates of (E}°°. Using a 4-term se- 
ries for the Aio basis set for this asymptotic series gave 
(E)°° = -2.879028802777 hartree which is in error by 
3.5 x 10~ 7 hartree. 

The ability to accurately predict (E)°° using eq. i|18|) 
has been tested for other values of A. Making the choice 
A = A 20 gave (E)°° = -2.87902876731919 hartree when 
the 6-term series was used to make the extrapolation. In 
summary, there is strong numerical evidence that eq. i|18|) 
correctly describes the convergence of the energy with N. 

The slower convergence of the electron-electron 5- 
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FIG. 6: The estimated exponents p$ as a function of 
for the LTO calculations of the He ground state (8). The 
different curves were obtained from the LTO basis sets with 
the exponents listed in Table U] 



FIG. 7: The extrapolated N —* oo limit for the He ground 
state energy obtained from the different asymptotic expan- 
sions. The energy sequence from the Aio was used. The exact 
s-wave energy as given by the calculation of Goldman et al 
4] is subtracted from the energy. 



function as AS N w A/TV 5 / 2 means that the ability to 
extrapolate to the N — > oo limit is even more impor- 
tant in obtaining accurate expectation values. Figure [S| 
shows the (S)°° estimates for the Aio basis while using 
eq. (f 1 3|) and eq. (|19|) to describe the large N limiting be- 
havior. It is noticed that the convergence improved as N p 
increased as long as N was sufficiently large. Choosing 
N > (M + 10) would seem to be sufficient for 2-term or 
3-term fits to eq. I|19|) . The specific numerical estimates 
of (S)°° for various extrapolations at N = 60 are given in 
Table UJ The 5-term fit gave (5) 00 = 0.1557637174 ajj. 
Given that {5) 6 ° for the A 10 basis was 0.155 772 7974 ag, 
the improvement in precision from the 5-term expansion 
corresponds to 4-5 orders of magnitude. This result is 
not specific to the Aio basis. Usage of the A20 basis re- 
sulted in an estimate of (5) 00 = 0.1557637175 ag. Given 
these results it would seem reasonable to assign a value 
of 0.1557637174(2) af, to (S). This is the value that was 
adopted as the "exact" value when plotting Figures|3|and 

El 

The more sophisticated extrapolations shown in Figure 
El exhibit even-odd fluctuations at the lower values of N 
and led us to omit consideration of a 6-term fit. These 
fluctuations are believed to arise from structural features 
of the helium ground state for reasons outlined in the 
next section. 



VI. CONVERGENCE OF A NATURAL 
ORBITAL BASIS 

The asymptotic form of the energy for a wave function 
written in terms of its natural orbital decomposition 1, 
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FIG. 8: The extrapolated N — > 00 limit for the He ground 
state {5) 00 for the different extrapolation methods applied to 
the Aio basis as described in the text. 



0,0] has also been re-examined for s-wave helium. First 
a very large calculation of the ground state wave function 
with a basis of 70 LTOs (A = 11.10) was performed. The 
one electron density matrix was then diagonalized and 
the resulting natural orbitals were used to define a new 
orbital basis ordered in terms of decreasing ground state 
occupancy. In its natural orbital form, the wave function 
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for a 1 S e state is written 

|*) =Y^ck Aij (ilHhHj\0O)<lH(Ti)Mra) ■ (20) 

i 

The natural orbital expansion is usually ordered in terms 
of decreasing \di\. 

Table ITTll gives (E) and (5) for the sequence of increas- 
ingly larger NO expansions. For these calculations, the 
generated NOs were added successively and (Et) and 
(6) N computed once the hamiltonian was diagonalized. 
The calculations were taken up to a maximum NO ex- 
pansion length of 20. The LTO basis of dimension 70 was 
not large enough to give a precise representation of the 
NOs beyond that point. The energies in the table are 
expected to be accurate estimates of the "exact" NO en- 
ergy for all digits with the possible exception of the last 
two. The energies in Table lllll are slightly lower than 
the previous tabulations of the s-wave NO energies by 
Carroll et al and Goldman 0. We treat the NO 
orbitals merely as a particularly optimal set of orbitals 
to input into a CI calculations. So unlike Carroll et al 
and Goldman the configuration space is not restricted to 
only include </>i(ri)^ i (r 2 ) type configurations. It should 
be noted that we have also done some calculations us- 
ing the pure NO configuration space and when this is 
done the energies agree with those of Carroll et al and 
Goldman to all digits. The (8) values in Table ITTT1 are 
expected to approximate those of the "exact" basis to 
about to 10 digits. 

Figure [5] shows the variation of pe and ps versus 1 /N 
for a sequence of increasingly larger NO calculations up 
to N = 20. The visual inspection of the px vs 1/N 
curve immediately suggests that pe — 6 + A/N + . . . and 
ps = 4 + A/N + . . .. The supposition has been confirmed 
by doing fits to the asymptotic form 



A' „ 



. \ - Pi 



(21) 



for increasingly larger values of N p . The results for po and 
pi are given in Table ITT1 The present calculations give po 
values of 5.992, 5.997 and 5.996 (the 5-term series which 
gave po = 5.967 is likely to be more susceptible to small 
imperfections in the NOs). A least-squares fit to the 
function pe = Po +pi/N t over the N £ [11,20] interval 
gavepo = 6.0005 and t = 1.070. A least-squares fit to the 
function ps = po + P\/N t over the N £ [11,20] interval 
gave po — 3.992 and t = 0.9174. A small even-odd ripple 
was present in the ps vs N graph. 

The linear variation of px vs 1 /N indicates the asymp- 
totic series 



AE N 
A5 N 



Ae Be Ce 

iV6 7V« 

As_ Bs_ Cs_ 

N 4 N 5 N 6 



(22) 
(23) 




for the variation of the AE N and A6 N with N. The 
0(N~ 6 ) variation of AE N has been known since the 



FIG. 9: The estimated exponents pe, Ps and p r as a function 
of for the NO expansion of He ground state. The value of 
p r shows strong even-odd oscillations. 



work of Carroll et al 4]. The present calculations give a 
more precise determination of the exponent. The value 
of p ps 5.7 previously reported by Goldman [l3T | can be 
discounted. Besides giving the leading order term with 
increased precision, the present calculations also demon- 
strate the order of the next term in the asymptotic se- 
ries is 0(N~ 7 ) and thereby strengthen the justification 
of asymptotic series based on expansions of the principal 
quantum number [T^ . 

The best NO estimate of (E)°° in Table |TT| was about 
an order of magnitude less accurate than that obtained 
from the LTO calculation. This was not unexpected since 
the values of (E) N in Table ITTTI are not the "exact" NO 
energies, merely very good estimates of these energies. 
Also, use of the NO basis inevitably means a more com- 
plicated calculation, and so it is more likely to be af- 
fected by round-off errors and discretization errors in the 
numerical quadratures. 

The coefficient of the leading order term for AE was 
-0.2999 hartree which is a bit larger in magnitude than 
the value initially given by Carroll et al [lj, namely - 
0.24 hartree. This level of agreement is acceptable given 
the fact that Carroll et al actually use a slightly differ- 
ent Ae/(N — 1/2) ~ 6 functional form (and do not allow 
for higher order terms) and extract the value of Ae at 
N w 10 which is too low to extract the asymptotic value 
of Ae (the value of Ae varies by more than 10% be- 
tween N = 10 and N = 20 for a single-term asymptotic 
formula). The precision of the Carroll et al calculation 
is also less than that of the present calculation (they ob- 
tained -2.879028765 hartree as their variational limit). 

The leading order term for the variation of AS N with 
N was 0(N~ 4 ). This dependence is consistent with ear- 
lier work of Halkier et al 20] who found that the variation 
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of A5 X with X to be 0(X~ 2 ) where X is the principal 
quantum number of the natural orbital. When analyz- 
ing this set of data it was discovered that there were 
regular fluctuations in the derived parameters as a func- 
tion of N as the analysis was made more sophisticated. 
When the 3-term approximation to eq. Ijl7|l was used 
the value of po oscillated between 3.94 and 4.06 depend- 
ing on whether N was even or odd (the oscillations in 
Pi were more marked). The oscillations became larger 
for the 4-term fit, here it was found that the po typ- 
ically flipped between 3.4 and 4.6. The actual values 
given in Table UTI were obtained by weighted average, e.g. 
po = 0.25p a (N = 19) + 0.50p (^ = 20) + 0.25p (^ = 21). 

These oscillations are most likely due to the physical 
properties of the basis set expansion of the He ground 
state. It has been known for a long time that treating 
the two electrons as an inner and outer electron can lead 
to a better description of the radial correlations @, • 
With the electrons having a tendency to separate into 
inner and outer radial orbitals the possibility does exist 
that achieving this separation might be slightly easier or 
harder if there are an even or odd number of NOs. It must 
be recalled that NOs themselves are objects that depend 
on the electron dynamics. The convergence of the mean 
electron-nucleus distance, i.e. (r) was also examined and 
the convergence pattern was quite irregular. Defining p r 
using eq. lfT4"jl results in the strongly oscillating plot of 
p r vs 1/N observed in Figure |5J The oscillations disap- 
pear if (r) for only even N (or only odd N) are used 
in a slightly modified version of eq. Ijl4|l and one finds 
the leading order term in Ar N is 0(N~ 6 ). It should be 
noted that similar even-odd fluctuations have also been 
observed in h igh precision calculations using correlated 
basis sets ^2 HE H3 ■ 

The asymptotic behavior of the natural orbital configu- 
ration coefficients were also determined. The coefficients 
are the di in eq. I|20l) . Assuming that the di scale as an 
inverse power series, dj w A^/i^Q gives 



In 



In 



i - 1 



(24) 



A fit of p to i using the formula 



PNO = PO + — + -7 , 



(25) 



gave values of po that ranged from 3.998 to 4.003 for suc- 
cessive fits to the 3 previous values for i-values between 
12 and 20 for the A60 basis. It was found that 



0.362 0.589 1.492 
—1— H nr- H tt— 



(26) 



at % = 20. Carroll et al obtained the result di rs ^rfjgp 



TABLE III: The term by term energy and (S) -function (in 
expectation values for the NO basis. 



N 






(E) 


N 








(S) 






1 


-2 


.861 


531 


101 


7265 


0. 


190 


249 


652 


529 


2 


-2 


.877 


929 


200 


9378 


0. 


161 


369 


548 


453 


3 


-2 


.878 


844 


196 


5241 


0. 


157 


747 


352 


993 


4 


-2 


.878 


980 


288 


7909 


0. 


156 


661 


432 


897 


5 


-2 


.879 


012 


046 


8823 


0. 


156 


240 


259 


959 


6 


-2 


.879 


021 


844 


0177 


0. 


156 


045 


264 


861 


7 


-2 


.879 


025 


501 


6647 


0. 


155 


943 


428 


536 


8 


-2 


.879 


027 


069 


6581 


0. 


155 


885 


238 


797 


9 


-2 


.879 


027 


815 


8906 


0. 


155 


849 


654 


582 


10 


-2 


.879 


028 


201 


2549 


0. 


155 


826 


694 


126 


11 


-2 


.879 


028 


413 


7401 


0. 


155 


811 


228 


411 


12 


-2 


.879 


028 


537 


3691 


0. 


155 


800 


434 


791 


13 


-2 


.879 


028 


612 


5987 


0. 


155 


792 


675 


857 


14 


-2 


.879 


028 


660 


1493 


0. 


155 


786 


956 


341 


15 


-2 


.879 


028 


691 


2007 


0. 


155 


782 


648 


350 


16 


-2 


.879 


028 


712 


0592 


0. 


155 


779 


342 


132 


17 


-2 


.879 


028 


726 


4219 


0. 


155 


776 


762 


792 


18 


-2 


.879 


028 


736 


5303 


0. 


155 


774 


721 


119 


19 


-2 


.879 


028 


743 


7842 


0. 


155 


773 


084 


071 


20 


-2 


.879 


028 


749 


0810 


0. 


155 


771 


756 


198 



VII. CONVERGENCE OF AN OPTIMIZED 
BASIS 

In this section the convergence properties of the LTO 
basis which is energy optimized at each N are studied. 
Developing the sequence of exponents Am that gave the 
lowest energy for a LTO basis of dimension M was te- 
dious. Defining 6(E) and S(S) as the differences in (E) 
and (5) arising from an imprecisely known Amj one has 
the relations 



5(E) 
5(5) 



A{5\f 
B(5\) . 



(27) 
(28) 



1] 



The quadratic dependence of 5(E) with respect to SA 
does make it easier to generate the sequence of (E) N val- 
ues. But this quadratic dependence upon 5X does make 
it harder to determine Am since the energy only depends 
weakly on A in the vicinity of the minimum. Since 5(5) 
depends linearly on 5A, any imprecision in Am impacts 
the precision of the (5) N sequence more severely. 

Some specific data can be used to put this in perspec- 
tive. The Am for M = 1, . . . , 30 have been determined 
to a precision for 10~ 6 for the calculations reported in 
this section. These gave an energy that was accurate to 
10 -18 hartree for the M = 15 calculation, but (5) was 
only known to a precision of 10 -11 ajj. Determination of 
(5) to a precision of 10 -15 Qq would require fixing Am 
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with an accuracy of 10~ 10 which would necessitate an 
energy given to a accuracy of 10~ 26 hartree. 

The behavior of pe and ps vs N was sufficiently com- 
plicated that an initial least squares fit to the equation 
p = Po + pi/N* was performed for N G [18,30]. The 
results of the fit gave 



Pe 
Ps 



= 5.6562 



= 3.8093 



15.69 

^2.7326 

0.3832 

^yO.5438 



(29) 
(30) 



The distinctive aspect about the fit is the difference in the 
leading terms of the inverse power series for pe and ps- 
Figure ITUl shows that variation of pe for the optimized 
LTO basis as a function of 1/N 2 - 7326 up to N = 30. The 
plot of ps is tending to curl up for the smallest values of 
jyjy-2.7326 b ecause [ t j s not linear in 1/N 2 - 7326 . 




0.001 0.002 

1/N 27326 



FIG. 10: The estimated exponents pe and pg as a function 
of JV -2 ' 7326 for the optimized basis. The variation of pg with 
TV -2,7326 is not expected to be linear. 

Another notable feature of Figure EH were the oscilla- 
tions in pe and ps for even and odd values of N. Oscil- 
lations in ps were previously seen for the NO sequence 
but the ps oscillations in Figure ITUl are more pronounced 
than those in Figure El Some of the values in Table |Tj] 
were given using the 3-point averaging used previously 
for the Natural Orbital sequence. 

The asymptotic analysis to determine the variational 
limits were performed with the following series 



AE N 
AS N 



A, 



Br 



jy5.6562 
As 



^y8.3888 

Bs 



^3.8093 1 ^4.3531 



(31) 
(32) 



The results of the analysis are given in Table [TJ] The 
energy is predicted with an accuracy of 10~ 10 hartree 



while ($) 00 is given to an accuracy of 10 -8 <Zq. Equations 
(|31|l and (|32|l were not worth extending to include more 
terms. The power of the next term in eq. H31JI is not 
obvious (refer to the Appendix A) and the oscillations 
in ps to a certain extent negate the value of extending 
eq. (|32|) to include additional terms (even if we knew 
what those terms were!). 



VIII. SUMMARY AND CONCLUSIONS 

Results of a sequence of CI calculations of the He 
ground state with an t = basis have been presented. 
This can be regarded as the simplest model of a real 
atom that has a correlation cusp. The energy depen- 
dence of the LTO basis was AE w 0{N~ 7 / 2 ). This 
rather slow convergence rate can be improved by fit- 
ting a succession of (E) values to the inverse power 
series AE J = Ae/N^ 7 / 2 + B E /N~ 8 / 2 + ... and esti- 
mating the N — > oo limit. It ultimately proved possible 
after adopting quadruple precision arithmetic, to repro- 
duce the known energy in this model to an accuracy of 
1 x 10 -14 hartree. The specific choice of the asymptotic 
series should be regarded as conjecture supported by nu- 
merical evidence. More definite proof would require the 
calculations to be extended to N > 100. The common 
exponent of the LTO basis should not be chosen to op- 
timize the energy for the largest calculation since this 
results in a distorted convergence pattern. In effect, op- 
timizing the LTO exponent for N LTOs, and then using 
the (A-3), (N-2), (N-l) and N energies to determine the 
coefficients of a 3-term expansion to eq. (|18fl will give an 
inaccurate estimate of the energy correction needed to 
achieve the variational limit. Any extrapolation would 
seem to require that N (the number of LTOs) should ex- 
ceed M (the basis dimension at which A was optimized) 
by about ten or more. This conclusion holds for both the 
energy and electron-electron 5-function. The very slow 
0(1/N 5 / 2 ) convergence of (S) was also circumvented by 
the use of the N — > oo corrections. 

The examinations of the convergence rate for an NO 
basis set revealed a faster convergence. The NO basis 
converged as 0(N~ 6 ) with the next term being 0(N~ 7 ). 
The present determinations of the convergence rates are 
more rigorous than those of Carroll et al lj. One sur- 
prising result was the slight even-odd oscillation in the 
convergence of the inter-electronic 5-function. Examina- 
tion of the (r) revealed noticeable even-odd oscillations 
in p r . The presence of these ripples could complicate de- 
termination of the variational limit of expectation values 
other than the energy. It was possible to extrapolate the 
energy of a 20 orbital NO basis to the variational limit 
with an accuracy of about 10~ 12 hartree. 

The convergence rate of the optimized LTO basis was 
0(N- 5 - 6562 ) with the next term being 0(A~ 8 - 3888 ). The 
degree of uncertainty in both of these exponents is much 
larger than for the fixed A LTO sequence or the NO 
sequence. The extremely tedious nature of the A op- 
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timization, combined with the lack of knowledge about 
the nature of the asymptotic series beyond the first two 
terms, make this extrapolation a less attractive proposi- 
tion. The noticeable even-odd oscillation in ps and even 
Pe further render the method even more unattractive. 
The implications of this behavior are not confined to the 
present work. For example, it is likely that correlated 
exponential basis sets composed of functions with 



£(ri,r 2 ,r 12 ) 



^exp(-Ar 1 )exp(-Ar 2 ) (33) 



could also exhibit complicated convergence patterns since 
A is often energy optimized as the basis dimension is in- 
creased in size |2£}. Consequently, it would not be sur- 
prising for estimates of the N — > oo energy correction 
for variational calculations on systems using a Hylleraas 
basis to be unreliable. For example, Yan and cowork- 
ers have estimated the variational limit in a high preci- 
sion calculation of PsH using a Hylleraas type basis [3l| . 
Their estimated energy correction for the PsH ground 
state energy (only 9.6 x 10~ 8 hartree) was too small by 
at least a factor of three [s3 • 

One of the main motivations for the present study 
was to gain insight into solving the problems associ- 
ated with the very slow convergence of CI calculations 
for mixed electron-positron systems 0, 0, |2jl |3j| . In 
effect, the problem is to determine the complete basis 
set limit 0, 0, EH for these exotic systems. The slow 
0{N~ 7 / 2 ) convergence of the energy for an LTO basis 
set is greatly improved by the adoption of extrapolation 
schemes. Using the N = 10 energy for the Aio basis and 
the best extrapolation of the N = 60 calculation in Table 
ITTI as two reference points, one deduces an effective con- 
vergence rate of 0(N~ 10 ). The penalty associated with 
the use of the extrapolation formulae is the necessity to 
use quadruple precision arithmetic if 3 or more terms are 
retained in the inverse power series (note, a 3-term se- 
ries for AS N was numerically stable in double precision 
arithmetic). The need to use the quadruple precision 
arithmetic is caused by the very small size of the AE N 
increments and the impact of round-off error on the fit to 
the inverse power series. One somewhat ironic feature is 
that it is necessary to use a basis that is not energy op- 
timized so that the extrapolation to the variational limit 
can be done reliably. 
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APPENDIX A: ANALYSIS OF THE 
p-DEPENDENCE 

Let us demonstrate that an asymptotic series 



AX N = 44- 



A 



B 



Ni Ni +t 



(with C = B/A) leads to p = q + F/N l when p is defined 
from successive AX N increments by 



V AX N J I "' \N - 1 
Substituting Al" and AI*" 1 from eq. (|A"T|l gives 



(A2) 



p = In 




(A3) 



The logarithm in the numerator can be split into two 
terms 



fAX N - 1 \ , ( N 
In I — . —tt I = q In 



V AX N 



N - 1 



In 



1 + 



c 



(N-iy 



1+-^- 



(A4) 



The first term conveniently cancels with the denomina- 
tor to give q. The argument of the second term can be 
expanded 



1- 



c 



(N-iy 



A ' 



1 



c_ tc \ (, _c_ c 2 

+ Y* + Y*+iJ ^ ~Y* + Y2* 
tC 



N t+i 

Using ln(l +x) « x leads to 



(A5) 



In 



'II c 
. 1 + & 



tC 
jVt+j 



(A6) 



The denominator is simplified using \a.(N/(N — 1)) 
ln(l + 1/(JV - 1)) w 1/JV to finally give 



p = q- 



tC_ 

N*- 



(A7) 



as required. If eq. I|A1J) has successive terms where the 
power increments by t — 1 or t = 1/2 indefinitely, then 
this leads to a corresponding series, eq. I|A7|) that also 
have powers that respectively increment by 1 or 1/2 in- 
definitely This is not necessarily true for arbitrary t in 
eq. (JJTJ. 
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APPENDIX B: SCALING OF THE 2-ELECTRON 
INTEGRALS 

The most time-consuming part of the calculation was 
the generation of the electron-electron and annihilation 
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matrix elements. However, the expense of this was 
greatly reduced by generating an initial set of integrals 
for a given A, and then using a scaling factor to generate 
the integral lists for other values of A. 
The basic integral that has to be done is 

R(n a ,n b ,n c ,n d ,X) = J J dn dr 2 N a (X)N b (X) 

x N c (X)N d (X)f a {Xn)f b (Xr 2 ) 
x V(r u r 2 )f c (X ri )f d (Xr 2 ) (Bl) 

All integrals can be defined in terms of 
R(n a ,n b ,n c ,n d , X = 1). Consider the integral (|B1|) 
and make the transformation Xr = u. Therefore 
r*i = ui/X and r 2 = u 2 /X. Similarly dr\ — dui/X and 
dr 2 = du 2 1 X and therefore 

R(n a: n b ,n c ,n d ,X) = J J du x du 2 N a {X)N b (X) 

x N c (X)N d {X)f a ( Ul )f b (u 2 ) 

x V( ri ,r 2 )f c ( Ul )f d (u 2 ) (B2) 



From eq. ©, N a (X) = X 1 / 2 N a (X = 1), so 

R(n a ,n b ,n c ,n d: A) = J J dux du 2 N a {l)N b {\) 

x N c (l)N d (l)f a ( Ul )f b (u 2 ) 

x V(ri,ra)/ C (ui)/,,(U2) (B3) 

The scaling for the electron-electron repulsion integral is 
|r a - r 2 | _1 = A|ui - u 2 | _1 . Hence 

R{n a , n b , n c , n d , X) = XR{n ai n bl n c , n d , 1) , (B4) 

for the electron-electron integral. When the operator is 
the J-function, one uses the result (5(r 1 —r 2 ) = A(5(u!— u 2 ) 
to give 

R(n a , n b , n c , n d , X) = XR(n a ,n b , n c , n d , 1) . (B5) 
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